function handle = arrowhead(x,y,clr,Asize,head)

%  arrowhead   Draw a arrow head 

%      This program was written to place arrowheads on trajectories 
%      calculated using one of MATLAB's ode solvers.  However, the program
%      does not require the solvers, it only requires two points as shown
%      in the figure below.

%                (x2,y2)              (x2,y2)              (x2,y2)
%      head=1       o         head=2     o         head=3     o                      
%                  /|\                  /|\                  /|\                      
%                 / | \                / | \                / | \                  
%                /  |  \              /  o  \              /  |  \             
%               /   |   \            / / | \ \            / _-o-_ \               
%              o----o----o          o/   |   \o          o    |    o           
%                   |                    |                    |                 
%                   |                    |                    |
%           (x1,y1) o            (x1,y1) o            (x1,y1) o
      
%      arrowhead(x,y,color,Asize,head)  

%      required arguments:  the points x=[x1 x2] and y=[y1 y2]
%      optional arguments:  
%         color: a string, same format used in the plot command
%         Asize = [Asize(1) Asize(2)]: this rescales the arrowhead
%             Asize(1) = the length in changed by this factor
%             Asize(2) = the width in changed by this factor
%         head: an integer, this specifies the shape of the arrow head
%             head = 1 straight back (the default)
%             head = 2 slanted back
%             head = 3 curved back

%	 Examples:
%	 [t,y] = ode45(@rhs,[0 100],[0 1]);
%	 hold on
%	 plot(y(:,1),y(:,2))
%	 i=40; ii=i+3;
%	 arrowhead([y(i,1) y(ii,1)],[y(i,2) y(ii,2)]);
%	 arrowhead([y(i,1) y(ii,1)],[y(i,2) y(ii,2)],'r');
%	 arrowhead([y(i,1) y(ii,1)],[y(i,2) y(ii,2)],'k',[2 1],2);
%	 arrowhead([y(i,1) y(ii,1)],[y(i,2) y(ii,2)],[],[0.4 1.2]);
%	 arrowhead([y(i,1) y(ii,1)],[y(i,2) y(ii,2)],[],[],3);

%	 As indicated in the above three figures, the arrow head is centered on the 
%	 line through (x1,y1), (x2,y2), with the front of the arrow head at (x2,y2).  
%	 In the above examples these points are from the computed solution of an ODE.  
%	 In this context, (x1,y1) is the solution at time t1, while (x2,y2) is the 
%	 solution at time t2, where t1 < t2.  Depending on how picky you are about 
%	 the arrow head being centered on the solution curve, you can adjust t1 to 
%	 rotate the arrow's direction.

%	 The default length and width of the arrow head are determined by the height 
%	 and width of the currently active plot, or subplot, window.  These values can
%	 be adjusted using the Asize argument.    

handle = [];

% errors
if nargin < 2
	error('You need to specify two points');
end
if (length(x) ~= 2) || (length(y) ~= 2),
	error('x and y vectors must each have two components');
end
if (x(1) == x(2)) && (y(1) == y(2)),
	error('The two points are equal - cannot determine direction!');
end

% set the defaults
if nargin <= 2
	clr = 'b';
end
if nargin <= 3
	Asize = [1 1];
end
if nargin <= 4
	head = 1;
end

% check if variables left empty
if isempty(clr)
	clr = 'b'; 
end;
if isempty(Asize)
	Asize = [1 1];
end

% determine and remember the hold status, toggle if necessary
if ishold,
	WasHold = 1;
else
	WasHold = 0;
	hold on;
end

x1 = x(1); x2 = x(2); y1 = y(1); y2 = y(2);

% get the axis dimensions
Scales = get(gca,'Position');

% get ranges for axis 
TheAxis = axis;
Xlimit = abs(TheAxis(2) - TheAxis(1));
Ylimit = abs(TheAxis(4) - TheAxis(3));

% 2*Whead = width of arrow head
Whead = Asize(2)*min(Scales(3),Scales(4))/50;

% Lhead = length of arrow head
Lhead = 2.8*Whead*Asize(1);

% calculate the two tips of the arrow head
if x1==x2
	xa=x1+Xlimit*Whead/Scales(3);
	ya=y2-sign(y2-y1)*Ylimit/Scales(4);
	xb=x1-Xlimit*Whead/Scales(3);
	yb=ya;
else
	% calculate slope of arrow head
	s0=Scales(4)*Xlimit/(Scales(3)*Ylimit);
	m = s0*(y2 - y1)/(x2 - x1);

	% calculate center-point of wide end
	if x1>x2
		x3=x2+Xlimit*Lhead/(Scales(3)*sqrt(1+m^2));
	else
		x3=x2-Xlimit*Lhead/(Scales(3)*sqrt(1+m^2));
	end
	y3=y2+m*(x3-x2)/s0;

	% calculate the two tips of the arrow head
	xa=x3+Xlimit*m*Whead/(Scales(3)*sqrt(1+m^2));
	ya=y3-(xa-x3)/(s0*m);
	xb=x3-Xlimit*m*Whead/(Scales(3)*sqrt(1+m^2));
	yb=y3-(xb-x3)/(s0*m);
end;

% the polygon forming the arrow head
if head==1
	xd = [x2 xa xb];
	yd = [y2 ya yb];
elseif head==2
	LL=0.7*Lhead;
	if x1==x2
		x4=x2;
		y4=y2-sign(y2-y1)*Ylimit*LL/Scales(4);
	elseif x1>x2
		x4=x2+Xlimit*LL/(Scales(3)*sqrt(1+m^2));
		y4=y2+m*(x4-x2)/s0;
	else
		x4=x2-Xlimit*LL/(Scales(3)*sqrt(1+m^2));
		y4=y2+m*(x4-x2)/s0;
	end
	xd = [x2 xa x4 xb];
	yd = [y2 ya y4 yb];
else
	z2 = [x2 y2];
	z3 = [x3 y3];
	zb = [xb yb];
	za = [xa ya];

	z23=norm(z2-z3);
	zb3=norm(zb-z3);
	za3=norm(za-z3);
	z2b=norm(z2-zb);
	z2a=norm(z2-za);

	db=(zb-z3);
	da=(za-z3);
	d=(z2-z3)/z23;
	q=(zb-z2)/z2b;
	qq=(za-z2)/z2a;

	zi = z2 + 0.7*(z2-z3);
	z2i=norm(z2-zi);

	b=z23/z2i;
	a=z2i/zb3^b;
	aa=z2i/za3^b;

	nn=6;
	for i=1:nn
		zl=z3+db*(i-1)/(nn-1);
		zll=z3+da*(i-1)/(nn-1);
		k=z2b*norm(zl-z3)/zb3;
		kk=z2a*norm(zll-z3)/za3;
		zt=z2+k*q;
		ztt=z2+kk*qq;
		g=norm(zt-zl)-a*norm(zb-zl)^b;
		gg=norm(ztt-zll)-aa*norm(za-zll)^b;
		zm=zl+g*d;
		zmm=zll+gg*d;
		if i==1
			xd = [zm(1)];  yd=[zm(2)];
		else
			xd = [zmm(1) xd zm(1)];  yd=[zmm(2) yd zm(2)];
		end;
	end;
	xd = [x2 xd];
	yd = [y2 yd];
end

% draw the arrowhead
fill(xd,yd,clr,'EdgeColor',clr);

% restore original axe ranges and hold status
axis(TheAxis);
if ~WasHold,
	hold off
end

